#######################################################
## Plot the conditional counterfactual distributions ##
#######################################################
pdf(paste("CondDFQF",str.sample,str.rho,str.cb,str.robust,str.boot,".pdf",sep=""));

##############################
### Distribution functions ###
##############################
## All Distributions ##
plot(range(ys), c(0,1), type="n",col=1, xlab="Log of hourly wage", lwd=2, ylab="Probability", 
     main=paste("Conditional DF ",main.str, sep=""), sub=" ");

lines(ys, F2.cond, lwd=1, col = 'firebrick'  );
lines(ys, F2.cond.c1rho, lwd=1, col = 'orange3');
lines(ys, F2.cond.c1rhopi, lwd=1, col = 'dark green'  );
lines(ys, F2.cond.c1rhopibeta, lwd=1, col = 'steelblue4'  );
lines(ys, F1.cond, lwd=1, col = 'dark grey'  );

legend("bottomright",legend=c(legend.counter2,
                              as.expression(bquote(rho~"="~.(legend.counter1))),
                              as.expression(bquote(rho*","*pi~"="~.(legend.counter1))),
                              as.expression(bquote(rho*","*pi*","*beta~"="~.(legend.counter1))),
                              legend.counter1),
       col= c("firebrick","orange3","dark green","steelblue4","dark grey"),lty=1);


## All Distributions with CB ##
plot(range(ys), c(0,1), type="n",col=1, xlab="Log of hourly wage", lwd=2, ylab="Probability", 
     main=paste("Conditional DF ",main.str, sep=""), sub=" ");

polygon(c(ys,rev(ys)),c(ubound.cond[,1],rev(lbound.cond[,1])), density=100, border=F, col='firebrick1', lty = 1, lwd = 1);
lines(ys, F2.cond, lwd=1, col = 'firebrick'  );

polygon(c(ys,rev(ys)),c(ubound.cond[,2],rev(lbound.cond[,2])), density=100, border=F, col='orange', lty = 1, lwd = 1);
lines(ys, F2.cond.c1rho, lwd=1, col = 'orange3');

polygon(c(ys,rev(ys)),c(ubound.cond[,3],rev(lbound.cond[,3])), density=100, border=F, col='light green', lty = 1, lwd = 1);
lines(ys, F2.cond.c1rhopi, lwd=1, col = 'dark green'  );

polygon(c(ys,rev(ys)),c(ubound.cond[,4],rev(lbound.cond[,4])), density=100, border=F, col='light blue', lty = 1, lwd = 1);
lines(ys, F2.cond.c1rhopibeta, lwd=1, col = 'steelblue4'  );

polygon(c(ys,rev(ys)),c(ubound.cond[,5],rev(lbound.cond[,5])), density=100, border=F, col='light grey', lty = 1, lwd = 1);
lines(ys, F1.cond, lwd=1, col = 'dark grey'  );

legend("bottomright",legend=c(legend.counter2,
                              as.expression(bquote(rho~"="~.(legend.counter1))),
                              as.expression(bquote(rho*","*pi~"="~.(legend.counter1))),
                              as.expression(bquote(rho*","*pi*","*beta~"="~.(legend.counter1))),
                              legend.counter1),
       col= c("firebrick","orange3","dark green","steelblue4","dark grey"),lty=1);

## 1 & 2 ##
plot(range(ys), c(0,1), type="n",col=1, xlab="Log of hourly wage", lwd=2, ylab="Probability", 
     main=paste("Conditional DF ",main.str, sep=""), sub=" ");

polygon(c(ys,rev(ys)),c(ubound.cond[,1],rev(lbound.cond[,1])), density=100, border=F, col='firebrick1', lty = 1, lwd = 1);
lines(ys, F2.cond, lwd=1, col = 'firebrick'  );

polygon(c(ys,rev(ys)),c(ubound.cond[,2],rev(lbound.cond[,2])), density=100, border=F, col='orange', lty = 1, lwd = 1);
lines(ys, F2.cond.c1rho, lwd=1, col = 'orange3');


legend("bottomright",legend=c(legend.counter2,
                              as.expression(bquote(rho~"="~.(legend.counter1)))),
       col= c("firebrick","orange3"),lty=1);


## 2 & 3 ##
plot(range(ys), c(0,1), type="n",col=1, xlab="Log of hourly wage", lwd=2, ylab="Probability", 
     main=paste("Conditional DF ",main.str, sep=""), sub=" ");

polygon(c(ys,rev(ys)),c(ubound.cond[,2],rev(lbound.cond[,2])), density=100, border=F, col='orange', lty = 1, lwd = 1);
lines(ys, F2.cond.c1rho, lwd=1, col = 'orange3');

polygon(c(ys,rev(ys)),c(ubound.cond[,3],rev(lbound.cond[,3])), density=100, border=F, col='light green', lty = 1, lwd = 1);
lines(ys, F2.cond.c1rhopi, lwd=1, col = 'dark green'  );


legend("bottomright",legend=c(as.expression(bquote(rho~"="~.(legend.counter1))),
                              as.expression(bquote(rho*","*pi~"="~.(legend.counter1)))),
       col= c("orange3","dark green"),lty=1);


## 3 & 4 ##
plot(range(ys), c(0,1), type="n",col=1, xlab="Log of hourly wage", lwd=2, ylab="Probability", 
     main=paste("Conditional DF ",main.str, sep=""), sub=" ");

polygon(c(ys,rev(ys)),c(ubound.cond[,3],rev(lbound.cond[,3])), density=100, border=F, col='light green', lty = 1, lwd = 1);
lines(ys, F2.cond.c1rhopi, lwd=1, col = 'dark green'  );

polygon(c(ys,rev(ys)),c(ubound.cond[,4],rev(lbound.cond[,4])), density=100, border=F, col='light blue', lty = 1, lwd = 1);
lines(ys, F2.cond.c1rhopibeta, lwd=1, col = 'steelblue4'  );


legend("bottomright",legend=c(as.expression(bquote(rho*","*pi~"="~.(legend.counter1))),
                              as.expression(bquote(rho*","*pi*","*beta~"="~.(legend.counter1)))),
       col= c("dark green","steelblue4"),lty=1);


## 4 & 5 ##
plot(range(ys), c(0,1), type="n",col=1, xlab="Log of hourly wage", lwd=2, ylab="Probability", 
     main=paste("Conditional DF ",main.str, sep=""), sub=" ");

polygon(c(ys,rev(ys)),c(ubound.cond[,4],rev(lbound.cond[,4])), density=100, border=F, col='light blue', lty = 1, lwd = 1);
lines(ys, F2.cond.c1rhopibeta, lwd=1, col = 'steelblue4'  );

polygon(c(ys,rev(ys)),c(ubound.cond[,5],rev(lbound.cond[,5])), density=100, border=F, col='light grey', lty = 1, lwd = 1);
lines(ys, F1.cond, lwd=1, col = 'dark grey'  );

legend("bottomright",legend=c(as.expression(bquote(rho*","*pi*","*beta~"="~.(legend.counter1))),
                              legend.counter1),
       col= c("steelblue4","dark grey"),lty=1);


############################################
### Distribution functions with Joint CB ###
############################################
## All Distributions with CB ##
plot(range(ys), c(0,1), type="n",col=1, xlab="Log of hourly wage", lwd=2, ylab="Probability", 
     main=paste("Conditional DF with Joint CB ",main.str, sep=""), sub=" ");

polygon(c(ys,rev(ys)),c(ubound.cond.joint[,1],rev(lbound.cond.joint[,1])), density=100, border=F, col='firebrick1', lty=1, lwd=1);
lines(ys, F2.cond, lwd=1, col = 'firebrick'  );

polygon(c(ys,rev(ys)),c(ubound.cond.joint[,2],rev(lbound.cond.joint[,2])), density=100, border=F, col='orange', lty=1, lwd=1);
lines(ys, F2.cond.c1rho, lwd=1, col = 'orange3');

polygon(c(ys,rev(ys)),c(ubound.cond.joint[,3],rev(lbound.cond.joint[,3])), density=100, border=F, col='light green', lty=1, lwd=1);
lines(ys, F2.cond.c1rhopi, lwd=1, col = 'dark green'  );

polygon(c(ys,rev(ys)),c(ubound.cond.joint[,4],rev(lbound.cond.joint[,4])), density=100, border=F, col='light blue', lty=1, lwd=1);
lines(ys, F2.cond.c1rhopibeta, lwd=1, col = 'steelblue4'  );

polygon(c(ys,rev(ys)),c(ubound.cond.joint[,5],rev(lbound.cond.joint[,5])), density=100, border=F, col='light grey', lty=1, lwd=1);
lines(ys, F1.cond, lwd=1, col = 'dark grey'  );

legend("bottomright",legend=c(legend.counter2,
                              as.expression(bquote(rho~"="~.(legend.counter1))),
                              as.expression(bquote(rho*","*pi~"="~.(legend.counter1))),
                              as.expression(bquote(rho*","*pi*","*beta~"="~.(legend.counter1))),
                              legend.counter1),
       col= c("firebrick","orange3","dark green","steelblue4","dark grey"),lty=1);

## 1 & 2 ##
plot(range(ys), c(0,1), type="n",col=1, xlab="Log of hourly wage", lwd=2, ylab="Probability", 
     main=paste("Conditional DF with Joint CB ",main.str, sep=""), sub=" ");

polygon(c(ys,rev(ys)),c(ubound.cond.joint[,1],rev(lbound.cond.joint[,1])), density=100, border=F, col='firebrick1', lty = 1, lwd = 1);
lines(ys, F2.cond, lwd=1, col = 'firebrick'  );

polygon(c(ys,rev(ys)),c(ubound.cond.joint[,2],rev(lbound.cond.joint[,2])), density=100, border=F, col='orange', lty = 1, lwd = 1);
lines(ys, F2.cond.c1rho, lwd=1, col = 'orange3');


legend("bottomright",legend=c(legend.counter2,
                              as.expression(bquote(rho~"="~.(legend.counter1)))),
       col= c("firebrick","orange3"),lty=1);


## 2 & 3 ##
plot(range(ys), c(0,1), type="n",col=1, xlab="Log of hourly wage", lwd=2, ylab="Probability", 
     main=paste("Conditional DF with Joint CB ",main.str, sep=""), sub=" ");

polygon(c(ys,rev(ys)),c(ubound.cond.joint[,2],rev(lbound.cond.joint[,2])), density=100, border=F, col='orange', lty = 1, lwd = 1);
lines(ys, F2.cond.c1rho, lwd=1, col = 'orange3');

polygon(c(ys,rev(ys)),c(ubound.cond.joint[,3],rev(lbound.cond.joint[,3])), density=100, border=F, col='light green', lty = 1, lwd = 1);
lines(ys, F2.cond.c1rhopi, lwd=1, col = 'dark green'  );


legend("bottomright",legend=c(as.expression(bquote(rho~"="~.(legend.counter1))),
                              as.expression(bquote(rho*","*pi~"="~.(legend.counter1)))),
       col= c("orange3","dark green"),lty=1);


## 3 & 4 ##
plot(range(ys), c(0,1), type="n",col=1, xlab="Log of hourly wage", lwd=2, ylab="Probability", 
     main=paste("Conditional DF with Joint CB ",main.str, sep=""), sub=" ");

polygon(c(ys,rev(ys)),c(ubound.cond.joint[,3],rev(lbound.cond.joint[,3])), density=100, border=F, col='light green', lty = 1, lwd = 1);
lines(ys, F2.cond.c1rhopi, lwd=1, col = 'dark green'  );

polygon(c(ys,rev(ys)),c(ubound.cond.joint[,4],rev(lbound.cond.joint[,4])), density=100, border=F, col='light blue', lty = 1, lwd = 1);
lines(ys, F2.cond.c1rhopibeta, lwd=1, col = 'steelblue4'  );


legend("bottomright",legend=c(as.expression(bquote(rho*","*pi~"="~.(legend.counter1))),
                              as.expression(bquote(rho*","*pi*","*beta~"="~.(legend.counter1)))),
       col= c("dark green","steelblue4"),lty=1);


## 4 & 5 ##
plot(range(ys), c(0,1), type="n",col=1, xlab="Log of hourly wage", lwd=2, ylab="Probability", 
     main=paste("Conditional DF with Joint CB ",main.str, sep=""), sub=" ");

polygon(c(ys,rev(ys)),c(ubound.cond.joint[,4],rev(lbound.cond.joint[,4])), density=100, border=F, col='light blue', lty = 1, lwd = 1);
lines(ys, F2.cond.c1rhopibeta, lwd=1, col = 'steelblue4'  );

polygon(c(ys,rev(ys)),c(ubound.cond.joint[,5],rev(lbound.cond.joint[,5])), density=100, border=F, col='light grey', lty = 1, lwd = 1);
lines(ys, F1.cond, lwd=1, col = 'dark grey'  );

legend("bottomright",legend=c(as.expression(bquote(rho*","*pi*","*beta~"="~.(legend.counter1))),
                              legend.counter1),
       col= c("steelblue4","dark grey"),lty=1);




##########################
### Quantile Functions ###
##########################
plot(c(0,1), range(ys), type="n",col=1, xlab="Quantile Index", lwd=2, ylab="Log of hourly wage", 
     main=paste("Observed wage", main.str, sep = ""));

lines(F2.cond, ys, lwd=1, col = 'firebrick'  );
lines(F2.cond.c1rho, ys, lwd=1, col = 'orange3');
lines(F2.cond.c1rhopi, ys, lwd=1, col = 'dark green'  );
lines(F2.cond.c1rhopibeta, ys, lwd=1, col = 'steelblue4'  );
lines(F1.cond, ys, lwd=1, col = 'dark grey'  );

legend("bottomright",legend=c(legend.counter2,
                              as.expression(bquote(rho~"="~.(legend.counter1))),
                              as.expression(bquote(rho*","*pi~"="~.(legend.counter1))),
                              as.expression(bquote(rho*","*pi*","*beta~"="~.(legend.counter1))),
                              legend.counter1),
       col= c("firebrick","orange3","dark green","steelblue4","dark grey"),lty=1);

## With CB ##
plot(c(0,1), range(ys), type="n",col=1, xlab="Quantile Index", lwd=2, ylab="Log of hourly wage", 
     main=paste("Observed wage", main.str, sep = ""));

polygon(c(ubound.cond[,1],rev(lbound.cond[,1])),c(ys,rev(ys)), density=100, border=F, col='firebrick1', lty = 1, lwd = 1);
lines(F2.cond, ys, lwd=1, col = 'firebrick'  );

polygon(c(ubound.cond[,2],rev(lbound.cond[,2])),c(ys,rev(ys)), density=100, border=F, col='orange', lty = 1, lwd = 1);
lines(F2.cond.c1rho, ys, lwd=1, col = 'orange3');

polygon(c(ubound.cond[,3],rev(lbound.cond[,3])),c(ys,rev(ys)), density=100, border=F, col='light green', lty = 1, lwd = 1);
lines(F2.cond.c1rhopi, ys, lwd=1, col = 'dark green'  );

polygon(c(ubound.cond[,4],rev(lbound.cond[,4])),c(ys,rev(ys)), density=100, border=F, col='light blue', lty = 1, lwd = 1);
lines(F2.cond.c1rhopibeta, ys, lwd=1, col = 'steelblue4'  );

polygon(c(ubound.cond[,5],rev(lbound.cond[,5])),c(ys,rev(ys)), density=100, border=F, col='light grey', lty = 1, lwd = 1);
lines(F1.cond, ys, lwd=1, col = 'dark grey'  );

legend("bottomright",legend=c(legend.counter2,
                              as.expression(bquote(rho~"="~.(legend.counter1))),
                              as.expression(bquote(rho*","*pi~"="~.(legend.counter1))),
                              as.expression(bquote(rho*","*pi*","*beta~"="~.(legend.counter1))),
                              legend.counter1),
       col= c("firebrick","orange3","dark green","steelblue4","dark grey"),lty=1);

## diff ##
plot(c(0,1), range(ys), type="n",col=1, xlab="Quantile Index", lwd=2, ylab="Log of hourly wage", 
     main=paste("Observed wage", main.str, sep = ""));

polygon(c(ubound.cond[,1],rev(lbound.cond[,1])), c(ys,rev(ys)), density=100, border=F, col='firebrick1', lty = 1, lwd = 1);
lines(F2.cond, ys, lwd=1, col = 'firebrick'  );

polygon(c(ubound.cond[,5],rev(lbound.cond[,5])), c(ys,rev(ys)), density=100, border=F, col='light grey', lty = 1, lwd = 1);
lines(F1.cond, ys, lwd=1, col = 'dark grey'  );

legend("bottomright", legend=c(legend.counter2,legend.counter1), col= c("firebrick","dark grey"), lty=1);


## 1 & 2 ##
plot(c(0,1), range(ys), type="n",col=1, xlab="Quantile Index", lwd=2, ylab="Log of hourly wage", 
     main=paste("Observed wage", main.str, sep = ""));

polygon(c(ubound.cond[,1],rev(lbound.cond[,1])),c(ys,rev(ys)), density=100, border=F, col='firebrick1', lty = 1, lwd = 1);
lines(F2.cond, ys, lwd=1, col = 'firebrick'  );

polygon(c(ubound.cond[,2],rev(lbound.cond[,2])),c(ys,rev(ys)), density=100, border=F, col='orange', lty = 1, lwd = 1);
lines(F2.cond.c1rho, ys, lwd=1, col = 'orange3');

legend("bottomright",legend=c(legend.counter2,
                              as.expression(bquote(rho~"="~.(legend.counter1)))), 
                              col= c("firebrick","orange3"),lty=1);

## 2 & 3 ##
plot(c(0,1), range(ys), type="n",col=1, xlab="Quantile Index", lwd=2, ylab="Log of hourly wage", 
     main=paste("Observed wage", main.str, sep = ""));

polygon(c(ubound.cond[,2],rev(lbound.cond[,2])),c(ys,rev(ys)), density=100, border=F, col='orange', lty = 1, lwd = 1);
lines(F2.cond.c1rho, ys, lwd=1, col = 'orange3');

polygon(c(ubound.cond[,3],rev(lbound.cond[,3])),c(ys,rev(ys)), density=100, border=F, col='light green', lty = 1, lwd = 1);
lines(F2.cond.c1rhopi, ys, lwd=1, col = 'dark green'  );

legend("bottomright",legend=c(as.expression(bquote(rho~"="~.(legend.counter1))),
                              as.expression(bquote(rho*","*pi~"="~.(legend.counter1)))),
       col= c("orange3","dark green"),lty=1);

## 3 & 4 ##
plot(c(0,1), range(ys), type="n",col=1, xlab="Quantile Index", lwd=2, ylab="Log of hourly wage", 
     main=paste("Observed wage", main.str, sep = ""));

polygon(c(ubound.cond[,3],rev(lbound.cond[,3])),c(ys,rev(ys)), density=100, border=F, col='light green', lty = 1, lwd = 1);
lines(F2.cond.c1rhopi, ys, lwd=1, col = 'dark green'  );

polygon(c(ubound.cond[,4],rev(lbound.cond[,4])),c(ys,rev(ys)), density=100, border=F, col='light blue', lty = 1, lwd = 1);
lines(F2.cond.c1rhopibeta, ys, lwd=1, col = 'steelblue4'  );

legend("bottomright",
       legend=c(as.expression(bquote(rho*","*pi~"="~.(legend.counter1))),
                as.expression(bquote(rho*","*pi*","*beta~"="~.(legend.counter1)))),
       col= c("dark green","steelblue4"),lty=1);

## 4 & 5 ##
plot(c(0,1), range(ys), type="n",col=1, xlab="Quantile Index", lwd=2, ylab="Log of hourly wage", 
     main=paste("Observed wage", main.str, sep = ""));

polygon(c(ubound.cond[,4],rev(lbound.cond[,4])),c(ys,rev(ys)), density=100, border=F, col='light blue', lty = 1, lwd = 1);
lines(F2.cond.c1rhopibeta, ys, lwd=1, col = 'steelblue4'  );

polygon(c(ubound.cond[,5],rev(lbound.cond[,5])),c(ys,rev(ys)), density=100, border=F, col='light grey', lty = 1, lwd = 1);
lines(F1.cond, ys, lwd=1, col = 'dark grey'  );

legend("bottomright",legend=c(as.expression(bquote(rho~","*pi~","*beta~"="~.(legend.counter1))),
                                         legend.counter1), 
       col= c("steelblue4","dark grey"),lty=1);


# With Average Effects #
## 1 & 2 ##
plot(c(0,1), range(ys), type="n",col=1, xlab="Quantile Index", lwd=2, ylab="Log of hourly wage", 
     main=paste("Observed wage", main.str, sep = ""));

polygon(c(ubound.cond[,1],rev(lbound.cond[,1])),c(ys,rev(ys)), density=100, border=F, col='firebrick1', lty = 1, lwd = 1);
lines(F2.cond, ys, lwd=1, col = 'firebrick'  );

polygon(c(ubound.cond[,2],rev(lbound.cond[,2])),c(ys,rev(ys)), density=100, border=F, col='orange', lty = 1, lwd = 1);
lines(F2.cond.c1rho, ys, lwd=1, col = 'orange3');

abline(h=e.avg.F2.cond, col='firebrick');
abline(h=ubound.avg.cond[1], col='firebrick', lty=2);
abline(h=lbound.avg.cond[1], col='firebrick', lty=2);

abline(h=e.avg.F2.cond.c1rho, col='orange3');
abline(h=ubound.avg.cond[2], col='orange3', lty=2);
abline(h=lbound.avg.cond[2], col='orange3', lty=2);

legend("bottomright",legend=c(legend.counter2,as.expression(bquote(rho~"="~.(legend.counter1)))),
       col= c("firebrick","orange3"),lty=1);

## 2 & 3 ##
plot(c(0,1), range(ys), type="n",col=1, xlab="Quantile Index", lwd=2, ylab="Log of hourly wage", 
     main=paste("Observed wage", main.str, sep = ""));

polygon(c(ubound.cond[,2],rev(lbound.cond[,2])),c(ys,rev(ys)), density=100, border=F, col='orange', lty = 1, lwd = 1);
lines(F2.cond.c1rho, ys, lwd=1, col = 'orange3');

polygon(c(ubound.cond[,3],rev(lbound.cond[,3])),c(ys,rev(ys)), density=100, border=F, col='light green', lty = 1, lwd = 1);
lines(F2.cond.c1rhopi, ys, lwd=1, col = 'dark green'  );

abline(h=e.avg.F2.cond.c1rho, col='orange3');
abline(h=ubound.avg.cond[2], col='orange3', lty=2);
abline(h=lbound.avg.cond[2], col='orange3', lty=2);

abline(h=e.avg.F2.cond.c1rhopi, col='dark green');
abline(h=ubound.avg.cond[3], col='dark green', lty=2);
abline(h=lbound.avg.cond[3], col='dark green', lty=2);

legend("bottomright",legend=c(as.expression(bquote(rho~"="~.(legend.counter1))),
                              as.expression(bquote(rho~","*pi~"="~.(legend.counter1)))),
       col= c("orange3","dark green"),lty=1);

## 3 & 4 ##
plot(c(0,1), range(ys), type="n",col=1, xlab="Quantile Index", lwd=2, ylab="Log of hourly wage", 
     main=paste("Observed wage", main.str, sep = ""));

polygon(c(ubound.cond[,3],rev(lbound.cond[,3])),c(ys,rev(ys)), density=100, border=F, col='light green', lty = 1, lwd = 1);
lines(F2.cond.c1rhopi, ys, lwd=1, col = 'dark green'  );

polygon(c(ubound.cond[,4],rev(lbound.cond[,4])),c(ys,rev(ys)), density=100, border=F, col='light blue', lty = 1, lwd = 1);
lines(F2.cond.c1rhopibeta, ys, lwd=1, col = 'steelblue4'  );

abline(h=e.avg.F2.cond.c1rhopi, col='dark green');
abline(h=ubound.avg.cond[3], col='dark green', lty=2);
abline(h=lbound.avg.cond[3], col='dark green', lty=2);

abline(h=e.avg.F2.cond.c1rhopibeta, col='steelblue4');
abline(h=ubound.avg.cond[4], col='steelblue4', lty=2);
abline(h=lbound.avg.cond[4], col='steelblue4', lty=2);

legend("bottomright",legend=c(as.expression(bquote(rho*","*pi~"="~.(legend.counter1))),
                              as.expression(bquote(rho*","*pi*","*beta~"="~.(legend.counter1)))),
       col= c("dark green","steelblue4"),lty=1);

## 4 & 5 ##
plot(c(0,1), range(ys), type="n",col=1, xlab="Quantile Index", lwd=2, ylab="Log of hourly wage", 
     main=paste("Observed wage", main.str, sep = ""));

polygon(c(ubound.cond[,4],rev(lbound.cond[,4])),c(ys,rev(ys)), density=100, border=F, col='light blue', lty = 1, lwd = 1);
lines(F2.cond.c1rhopibeta, ys, lwd=1, col = 'steelblue4'  );

polygon(c(ubound.cond[,5],rev(lbound.cond[,5])),c(ys,rev(ys)), density=100, border=F, col='light grey', lty = 1, lwd = 1);
lines(F1.cond, ys, lwd=1, col = 'dark grey'  );

abline(h=e.avg.F2.cond.c1rhopibeta, col='steelblue4');
abline(h=ubound.avg.cond[4], col='steelblue4', lty=2);
abline(h=lbound.avg.cond[4], col='steelblue4', lty=2);

abline(h=e.avg.F1.cond, col='dark grey');
abline(h=ubound.avg.cond[5], col='dark grey', lty=2);
abline(h=lbound.avg.cond[5], col='dark grey', lty=2);

legend("bottomright",legend=c(as.expression(bquote(rho*","*pi*","*beta~"="~.(legend.counter1))),
                              legend.counter1),
       col= c("steelblue4","dark grey"),lty=1);


# With Heckman Decomposition #
## 1 & 2 ##
plot(c(0,1), range(ys), type="n",col=1, xlab="Quantile", lwd=2, ylab="Log of hourly wage", 
     main=paste("Conditional QF and Heckman, ",main.str, sep=""), sub=" ");

polygon(c(ubound.cond[,1],rev(lbound.cond[,1])),c(ys,rev(ys)), density=100, border=F, col='firebrick1', lty = 1, lwd = 1);
lines(F2.cond, ys, lwd=1, col = 'firebrick'  );

polygon(c(ubound.cond[,2],rev(lbound.cond[,2])),c(ys,rev(ys)), density=100, border=F, col='orange', lty = 1, lwd = 1);
lines(F2.cond.c1rho, ys, lwd=1, col = 'orange3');

lines(F2.Heck.cond, ys, lwd=1, col='firebrick');
lines(F2.Heck.cond.c1rho, ys, lwd=1, col='orange3');

legend("bottomright",legend=c(legend.counter2,as.expression(bquote(rho*"="~.(legend.counter1)))), col= c("firebrick","orange3"),lty=1);

## 2 & 3 ##
plot(c(0,1), range(ys), type="n",col=1, xlab="Quantile", lwd=2, ylab="Log of hourly wage", 
     main=paste("Conditional QF and Heckman, ",main.str, sep=""), sub=" ");

polygon(c(ubound.cond[,2],rev(lbound.cond[,2])),c(ys,rev(ys)), density=100, border=F, col='orange', lty = 1, lwd = 1);
lines(F2.cond.c1rho, ys, lwd=1, col = 'orange3');

polygon(c(ubound.cond[,3],rev(lbound.cond[,3])),c(ys,rev(ys)), density=100, border=F, col='light green', lty = 1, lwd = 1);
lines(F2.cond.c1rhopi, ys, lwd=1, col = 'dark green'  );

lines(F2.Heck.cond.c1rho, ys, lwd=1, col='orange3');
lines(F2.Heck.cond.c1rhopi, ys, lwd=1, col='dark green');

legend("bottomright",legend=c(as.expression(bquote(rho*"="~.(legend.counter1))),
                              as.expression(bquote(rho*","*pi~"="~.(legend.counter1)))),
       col= c("orange3","dark green"),lty=1);

## 3 & 4 ##
plot(c(0,1), range(ys), type="n",col=1, xlab="Quantile", lwd=2, ylab="Log of hourly wage", 
     main=paste("Conditional QF and Heckman, ",main.str, sep=""), sub=" ");

polygon(c(ubound.cond[,3],rev(lbound.cond[,3])),c(ys,rev(ys)), density=100, border=F, col='light green', lty = 1, lwd = 1);
lines(F2.cond.c1rhopi, ys, lwd=1, col = 'dark green'  );

polygon(c(ubound.cond[,4],rev(lbound.cond[,4])),c(ys,rev(ys)), density=100, border=F, col='light blue', lty = 1, lwd = 1);
lines(F2.cond.c1rhopibeta, ys, lwd=1, col = 'steelblue4'  );

lines(F2.Heck.cond.c1rhopi, ys, lwd=1, col='dark green');
lines(F2.Heck.cond.c1rhopibeta, ys, lwd=1, col='steelblue4');

legend("bottomright",legend=c(as.expression(bquote(rho*","*pi~"="~.(legend.counter1))),
                              as.expression(bquote(rho*","*pi*","*beta~"="~.(legend.counter1)))),
       col= c("dark green","steelblue4"),lty=1);

## 4 & 5 ##
plot(c(0,1), range(ys), type="n",col=1, xlab="Quantile", lwd=2, ylab="Log of hourly wage", 
     main=paste("Conditional QF and Heckman, ",main.str, sep=""), sub=" ");

polygon(c(ubound.cond[,4],rev(lbound.cond[,4])),c(ys,rev(ys)), density=100, border=F, col='light blue', lty = 1, lwd = 1);
lines(F2.cond.c1rhopibeta, ys, lwd=1, col = 'steelblue4'  );

polygon(c(ubound.cond[,5],rev(lbound.cond[,5])),c(ys,rev(ys)), density=100, border=F, col='light grey', lty = 1, lwd = 1);
lines(F1.cond, ys, lwd=1, col = 'dark grey'  );

lines(F2.Heck.cond.c1rhopibeta, ys, lwd=1, col='steelblue4');
lines(F1.Heck.cond, ys, lwd=1, col='dark grey');

legend("bottomright",legend=c(as.expression(bquote(rho*","*pi*","*beta~"="~.(legend.counter1))),
                              legend.counter1),
       col= c("steelblue4","dark grey"),lty=1);


########################################
### Quantile Functions with joint CB ###
########################################
plot(c(0,1), range(ys), type="n",col=1, xlab="Quantile", lwd=2, ylab="Log of hourly wage", 
     main=paste("Conditional QF with Joint CB, ",main.str, sep=""), sub=" ");

polygon(c(ubound.cond.joint[,1],rev(lbound.cond.joint[,1])),c(ys,rev(ys)), density=100, border=F, col='firebrick1', lty = 1, lwd = 1);
lines(F2.cond, ys, lwd=1, col = 'firebrick'  );

polygon(c(ubound.cond.joint[,2],rev(lbound.cond.joint[,2])),c(ys,rev(ys)), density=100, border=F, col='orange', lty = 1, lwd = 1);
lines(F2.cond.c1rho, ys, lwd=1, col = 'orange3');

polygon(c(ubound.cond.joint[,3],rev(lbound.cond.joint[,3])),c(ys,rev(ys)), density=100, border=F, col='light green', lty = 1, lwd = 1);
lines(F2.cond.c1rhopi, ys, lwd=1, col = 'dark green'  );

polygon(c(ubound.cond.joint[,4],rev(lbound.cond.joint[,4])),c(ys,rev(ys)), density=100, border=F, col='light blue', lty = 1, lwd = 1);
lines(F2.cond.c1rhopibeta, ys, lwd=1, col = 'steelblue4'  );

polygon(c(ubound.cond.joint[,5],rev(lbound.cond.joint[,5])),c(ys,rev(ys)), density=100, border=F, col='light grey', lty = 1, lwd = 1);
lines(F1.cond, ys, lwd=1, col = 'dark grey'  );

legend("bottomright",legend=c(legend.counter2,
                              as.expression(bquote(rho~"="~.(legend.counter1))),
                              as.expression(bquote(rho*","*pi~"="~.(legend.counter1))),
                              as.expression(bquote(rho*","*pi*","*beta~"="~.(legend.counter1))),
                              legend.counter1),
       col= c("firebrick","orange3","dark green","steelblue4","dark grey"),lty=1);

## diff ##
plot(c(0,1), range(ys), type="n",col=1, xlab="Quantile", lwd=2, ylab="Log of hourly wage", 
     main=paste("Conditional QF with Joint CB, ",main.str, sep=""), sub=" ");

polygon(c(ubound.cond.joint[,1],rev(lbound.cond.joint[,1])), c(ys,rev(ys)), density=100, border=F, col='firebrick1', lty = 1, lwd = 1);
lines(F2.cond, ys, lwd=1, col = 'firebrick'  );

polygon(c(ubound.cond.joint[,5],rev(lbound.cond.joint[,5])), c(ys,rev(ys)), density=100, border=F, col='light grey', lty = 1, lwd = 1);
lines(F1.cond, ys, lwd=1, col = 'dark grey'  );

legend("bottomright", legend=c(legend.counter2,legend.counter1), col= c("firebrick","dark grey"), lty=1);


## 1 & 2 ##
plot(c(0,1), range(ys), type="n",col=1, xlab="Quantile", lwd=2, ylab="Log of hourly wage", 
     main="Selection Sorting", sub=" ");

polygon(c(ubound.cond.joint[,1],rev(lbound.cond.joint[,1])),c(ys,rev(ys)), density=100, border=F, col='firebrick1', lty = 1, lwd = 1);
lines(F2.cond, ys, lwd=1, col = 'firebrick'  );

polygon(c(ubound.cond.joint[,2],rev(lbound.cond.joint[,2])),c(ys,rev(ys)), density=100, border=F, col='orange', lty = 1, lwd = 1);
lines(F2.cond.c1rho, ys, lwd=1, col = 'orange3');

legend("bottomright",legend=c(legend.counter2,
                              as.expression(bquote(rho~"="~.(legend.counter1)))), 
       col= c("firebrick","orange3"),lty=1);

## 2 & 3 ##
plot(c(0,1), range(ys), type="n",col=1, xlab="Quantile", lwd=2, ylab="Log of hourly wage", 
     main="Selection Structure", sub=" ");

polygon(c(ubound.cond.joint[,2],rev(lbound.cond.joint[,2])),c(ys,rev(ys)), density=100, border=F, col='orange', lty = 1, lwd = 1);
lines(F2.cond.c1rho, ys, lwd=1, col = 'orange3');

polygon(c(ubound.cond.joint[,3],rev(lbound.cond.joint[,3])),c(ys,rev(ys)), density=100, border=F, col='light green', lty = 1, lwd = 1);
lines(F2.cond.c1rhopi, ys, lwd=1, col = 'dark green'  );

legend("bottomright",legend=c(as.expression(bquote(rho~"="~.(legend.counter1))),
                              as.expression(bquote(rho*","*pi~"="~.(legend.counter1)))),
       col= c("orange3","dark green"),lty=1);

## 3 & 4 ##
plot(c(0,1), range(ys), type="n",col=1, xlab="Quantile", lwd=2, ylab="Log of hourly wage", 
     main="Wage Structure", sub=" ");

polygon(c(ubound.cond.joint[,3],rev(lbound.cond.joint[,3])),c(ys,rev(ys)), density=100, border=F, col='light green', lty = 1, lwd = 1);
lines(F2.cond.c1rhopi, ys, lwd=1, col = 'dark green'  );

polygon(c(ubound.cond.joint[,4],rev(lbound.cond.joint[,4])),c(ys,rev(ys)), density=100, border=F, col='light blue', lty = 1, lwd = 1);
lines(F2.cond.c1rhopibeta, ys, lwd=1, col = 'steelblue4'  );

legend("bottomright",legend=c(as.expression(bquote(rho*","*pi~"="~.(legend.counter1))),
                              as.expression(bquote(rho*","*pi*","*beta~"="~.(legend.counter1)))),
       col= c("dark green","steelblue4"),lty=1);

## 4 & 5 ##
plot(c(0,1), range(ys), type="n",col=1, xlab="Quantile", lwd=2, ylab="Log of hourly wage", 
     main="Composition", sub=" ");

polygon(c(ubound.cond.joint[,4],rev(lbound.cond.joint[,4])),c(ys,rev(ys)), density=100, border=F, col='light blue', lty = 1, lwd = 1);
lines(F2.cond.c1rhopibeta, ys, lwd=1, col = 'steelblue4'  );

polygon(c(ubound.cond.joint[,5],rev(lbound.cond.joint[,5])),c(ys,rev(ys)), density=100, border=F, col='light grey', lty = 1, lwd = 1);
lines(F1.cond, ys, lwd=1, col = 'dark grey'  );

legend("bottomright",legend=c(as.expression(bquote(rho*","*pi*","*beta~"="~.(legend.counter1))),
                              legend.counter1), 
       col= c("steelblue4","dark grey"),lty=1);


# With Average Effects #
## 1 & 2 ##
plot(c(0,1), range(ys), type="n",col=1, xlab="Quantile", lwd=2, ylab="Log of hourly wage", 
     main="Selection Sorting", sub=" ");

polygon(c(ubound.cond.joint[,1],rev(lbound.cond.joint[,1])),c(ys,rev(ys)), density=100, border=F, col='firebrick1', lty = 1, lwd = 1);
lines(F2.cond, ys, lwd=1, col = 'firebrick'  );

polygon(c(ubound.cond.joint[,2],rev(lbound.cond.joint[,2])),c(ys,rev(ys)), density=100, border=F, col='orange', lty = 1, lwd = 1);
lines(F2.cond.c1rho, ys, lwd=1, col = 'orange3');

abline(h=e.avg.F2.cond, col='firebrick');
abline(h=ubound.avg.cond.joint[1], col='firebrick', lty=2);
abline(h=lbound.avg.cond.joint[1], col='firebrick', lty=2);

abline(h=e.avg.F2.cond.c1rho, col='orange3');
abline(h=ubound.avg.cond.joint[2], col='orange3', lty=2);
abline(h=lbound.avg.cond.joint[2], col='orange3', lty=2);

legend("bottomright",legend=c(legend.counter2,
                              as.expression(bquote(rho~"="~.(legend.counter1)))),
       col= c("firebrick","orange3"),lty=1);

## 2 & 3 ##
plot(c(0,1), range(ys), type="n",col=1, xlab="Quantile", lwd=2, ylab="Log of hourly wage", 
     main="Selection Structure", sub=" ");

polygon(c(ubound.cond.joint[,2],rev(lbound.cond.joint[,2])),c(ys,rev(ys)), density=100, border=F, col='orange', lty = 1, lwd = 1);
lines(F2.cond.c1rho, ys, lwd=1, col = 'orange3');

polygon(c(ubound.cond.joint[,3],rev(lbound.cond.joint[,3])),c(ys,rev(ys)), density=100, border=F, col='light green', lty = 1, lwd = 1);
lines(F2.cond.c1rhopi, ys, lwd=1, col = 'dark green'  );

abline(h=e.avg.F2.cond.c1rho, col='orange3');
abline(h=ubound.avg.cond.joint[2], col='orange3', lty=2);
abline(h=lbound.avg.cond.joint[2], col='orange3', lty=2);

abline(h=e.avg.F2.cond.c1rhopi, col='dark green');
abline(h=ubound.avg.cond.joint[3], col='dark green', lty=2);
abline(h=lbound.avg.cond.joint[3], col='dark green', lty=2);

legend("bottomright",legend=c(as.expression(bquote(rho~"="~.(legend.counter1))),
                              as.expression(bquote(rho*","*pi~"="~.(legend.counter1)))),
       col= c("orange3","dark green"),lty=1);

## 3 & 4 ##
plot(c(0,1), range(ys), type="n",col=1, xlab="Quantile", lwd=2, ylab="Log of hourly wage", 
     main="Wage Structure", sub=" ");

polygon(c(ubound.cond.joint[,3],rev(lbound.cond.joint[,3])),c(ys,rev(ys)), density=100, border=F, col='light green', lty = 1, lwd = 1);
lines(F2.cond.c1rhopi, ys, lwd=1, col = 'dark green'  );

polygon(c(ubound.cond.joint[,4],rev(lbound.cond.joint[,4])),c(ys,rev(ys)), density=100, border=F, col='light blue', lty = 1, lwd = 1);
lines(F2.cond.c1rhopibeta, ys, lwd=1, col = 'steelblue4'  );

abline(h=e.avg.F2.cond.c1rhopi, col='dark green');
abline(h=ubound.avg.cond.joint[3], col='dark green', lty=2);
abline(h=lbound.avg.cond.joint[3], col='dark green', lty=2);

abline(h=e.avg.F2.cond.c1rhopibeta, col='steelblue4');
abline(h=ubound.avg.cond.joint[4], col='steelblue4', lty=2);
abline(h=lbound.avg.cond.joint[4], col='steelblue4', lty=2);

legend("bottomright",legend=c(as.expression(bquote(rho*","*pi~"="~.(legend.counter1))),
                              as.expression(bquote(rho*","*pi*","*beta~"="~.(legend.counter1)))),
       col= c("dark green","steelblue4"),lty=1);

## 4 & 5 ##
plot(c(0,1), range(ys), type="n",col=1, xlab="Quantile", lwd=2, ylab="Log of hourly wage", 
     main="Composition", sub=" ");

polygon(c(ubound.cond.joint[,4],rev(lbound.cond.joint[,4])),c(ys,rev(ys)), density=100, border=F, col='light blue', lty = 1, lwd = 1);
lines(F2.cond.c1rhopibeta, ys, lwd=1, col = 'steelblue4'  );

polygon(c(ubound.cond.joint[,5],rev(lbound.cond.joint[,5])),c(ys,rev(ys)), density=100, border=F, col='light grey', lty = 1, lwd = 1);
lines(F1.cond, ys, lwd=1, col = 'dark grey'  );

abline(h=e.avg.F2.cond.c1rhopibeta, col='steelblue4');
abline(h=ubound.avg.cond.joint[4], col='steelblue4', lty=2);
abline(h=lbound.avg.cond.joint[4], col='steelblue4', lty=2);

abline(h=e.avg.F1.cond, col='dark grey');
abline(h=ubound.avg.cond.joint[5], col='dark grey', lty=2);
abline(h=lbound.avg.cond.joint[5], col='dark grey', lty=2);

legend("bottomright",legend=c(as.expression(bquote(rho*","*pi*","*beta~"="~.(legend.counter1))),
                              legend.counter1),
       col= c("steelblue4","dark grey"),lty=1);



############################################
### Quantile Functions with Bootstrap CB ###
############################################
plot(c(0,1), range(ys), type="n",col=1, xlab="Quantile", lwd=2, ylab="Log of hourly wage", 
     main=paste("Conditional QF with Joint Bootstrap CB, ",main.str, sep=""), sub=" ");

polygon(c(seq(l.thres,u.thres,by.thres),seq(u.thres,l.thres,-by.thres)),
        c(ubound.QFcond.joint[,1],rev(lbound.QFcond.joint[,1])), density=100, border=F, col='firebrick1', lty = 1, lwd = 1);
lines(F2.cond, ys, lwd=1, col = 'firebrick'  );

polygon(c(seq(l.thres,u.thres,by.thres),seq(u.thres,l.thres,-by.thres)),
        c(ubound.QFcond.joint[,2],rev(lbound.QFcond.joint[,2])), density=100, border=F, col='orange', lty = 1, lwd = 1);
lines(F2.cond.c1rho, ys, lwd=1, col = 'orange3');

polygon(c(seq(l.thres,u.thres,by.thres),seq(u.thres,l.thres,-by.thres)),
        c(ubound.QFcond.joint[,3],rev(lbound.QFcond.joint[,3])), density=100, border=F, col='light green', lty = 1, lwd = 1);
lines(F2.cond.c1rhopi, ys, lwd=1, col = 'dark green'  );

polygon(c(seq(l.thres,u.thres,by.thres),seq(u.thres,l.thres,-by.thres)),
        c(ubound.QFcond.joint[,4],rev(lbound.QFcond.joint[,4])), density=100, border=F, col='light blue', lty = 1, lwd = 1);
lines(F2.cond.c1rhopibeta, ys, lwd=1, col = 'steelblue4'  );

polygon(c(seq(l.thres,u.thres,by.thres),seq(u.thres,l.thres,-by.thres)),
        c(ubound.QFcond.joint[,5],rev(lbound.QFcond.joint[,5])), density=100, border=F, col='light grey', lty = 1, lwd = 1);
lines(F1.cond, ys, lwd=1, col = 'dark grey'  );

legend("bottomright",legend=c(legend.counter2,
                              as.expression(bquote(rho~"="~.(legend.counter1))),
                              as.expression(bquote(rho*","*pi~"="~.(legend.counter1))),
                              as.expression(bquote(rho*","*pi*","*beta~"="~.(legend.counter1))),
                              legend.counter1),
       col= c("firebrick","orange3","dark green","steelblue4","dark grey"),lty=1);

## diff ##
plot(c(0,1), range(ys), type="n",col=1, xlab="Quantile", lwd=2, ylab="Log of hourly wage", 
     main=paste("Conditional QF with Joint Bootstrap CB, ",main.str, sep=""), sub=" ");

polygon(c(seq(l.thres,u.thres,by.thres),seq(u.thres,l.thres,-by.thres)),
        c(ubound.QFcond.joint[,1],rev(lbound.QFcond.joint[,1])), density=100, border=F, col='firebrick1', lty = 1, lwd = 1);
lines(F2.cond, ys, lwd=1, col = 'firebrick'  );

polygon(c(seq(l.thres,u.thres,by.thres),seq(u.thres,l.thres,-by.thres)),
        c(ubound.QFcond.joint[,5],rev(lbound.QFcond.joint[,5])), density=100, border=F, col='light grey', lty = 1, lwd = 1);
lines(F1.cond, ys, lwd=1, col = 'dark grey'  );

legend("bottomright", legend=c(legend.counter2,legend.counter1), col= c("firebrick","dark grey"), lty=1);


## 1 & 2 ##
plot(c(0,1), range(ys), type="n",col=1, xlab="Quantile", lwd=2, ylab="Log of hourly wage", 
     main=paste("Conditional QF with Joint Bootstrap CB, ",main.str, sep=""), sub=" ");

polygon(c(seq(l.thres,u.thres,by.thres),seq(u.thres,l.thres,-by.thres)),
        c(ubound.QFcond.joint[,1],rev(lbound.QFcond.joint[,1])), density=100, border=F, col='firebrick1', lty = 1, lwd = 1);
lines(F2.cond, ys, lwd=1, col = 'firebrick'  );

polygon(c(seq(l.thres,u.thres,by.thres),seq(u.thres,l.thres,-by.thres)),
        c(ubound.QFcond.joint[,2],rev(lbound.QFcond.joint[,2])), density=100, border=F, col='orange', lty = 1, lwd = 1);
lines(F2.cond.c1rho, ys, lwd=1, col = 'orange3');

legend("bottomright",legend=c(legend.counter2,
                              as.expression(bquote(rho~"="~.(legend.counter1)))), 
       col= c("firebrick","orange3"),lty=1);

## 2 & 3 ##
plot(c(0,1), range(ys), type="n",col=1, xlab="Quantile", lwd=2, ylab="Log of hourly wage", 
     main=paste("Conditional QF with Joint Bootstrap CB, ",main.str, sep=""), sub=" ");

polygon(c(seq(l.thres,u.thres,by.thres),seq(u.thres,l.thres,-by.thres)),
        c(ubound.QFcond.joint[,2],rev(lbound.QFcond.joint[,2])), density=100, border=F, col='orange', lty = 1, lwd = 1);
lines(F2.cond.c1rho, ys, lwd=1, col = 'orange3');

polygon(c(seq(l.thres,u.thres,by.thres),seq(u.thres,l.thres,-by.thres)),
        c(ubound.QFcond.joint[,3],rev(lbound.QFcond.joint[,3])), density=100, border=F, col='light green', lty = 1, lwd = 1);
lines(F2.cond.c1rhopi, ys, lwd=1, col = 'dark green'  );

legend("bottomright",legend=c(as.expression(bquote(rho~"="~.(legend.counter1))),
                              as.expression(bquote(rho*","*pi~"="~.(legend.counter1)))),
       col= c("orange3","dark green"),lty=1);

## 3 & 4 ##
plot(c(0,1), range(ys), type="n",col=1, xlab="Quantile", lwd=2, ylab="Log of hourly wage", 
     main=paste("Conditional QF with Joint Bootstrap CB, ",main.str, sep=""), sub=" ");

polygon(c(seq(l.thres,u.thres,by.thres),seq(u.thres,l.thres,-by.thres)),
        c(ubound.QFcond.joint[,3],rev(lbound.QFcond.joint[,3])), density=100, border=F, col='light green', lty = 1, lwd = 1);
lines(F2.cond.c1rhopi, ys, lwd=1, col = 'dark green'  );

polygon(c(seq(l.thres,u.thres,by.thres),seq(u.thres,l.thres,-by.thres)),
        c(ubound.QFcond.joint[,4],rev(lbound.QFcond.joint[,4])), density=100, border=F, col='light blue', lty = 1, lwd = 1);
lines(F2.cond.c1rhopibeta, ys, lwd=1, col = 'steelblue4'  );

legend("bottomright",legend=c(as.expression(bquote(rho*","*pi~"="~.(legend.counter1))),
                              as.expression(bquote(rho*","*pi*","*beta~"="~.(legend.counter1)))),
       col= c("dark green","steelblue4"),lty=1);

## 4 & 5 ##
plot(c(0,1), range(ys), type="n",col=1, xlab="Quantile", lwd=2, ylab="Log of hourly wage", 
     main=paste("Conditional QF with Joint Bootstrap CB, ",main.str, sep=""), sub=" ");

polygon(c(seq(l.thres,u.thres,by.thres),seq(u.thres,l.thres,-by.thres)),
        c(ubound.QFcond.joint[,4],rev(lbound.QFcond.joint[,4])), density=100, border=F, col='light blue', lty = 1, lwd = 1);
lines(F2.cond.c1rhopibeta, ys, lwd=1, col = 'steelblue4'  );

polygon(c(seq(l.thres,u.thres,by.thres),seq(u.thres,l.thres,-by.thres)),
        c(ubound.QFcond.joint[,5],rev(lbound.QFcond.joint[,5])), density=100, border=F, col='light grey', lty = 1, lwd = 1);
lines(F1.cond, ys, lwd=1, col = 'dark grey'  );

legend("bottomright",legend=c(as.expression(bquote(rho*","*pi*","*beta~"="~.(legend.counter1))),
                              legend.counter1), 
       col= c("steelblue4","dark grey"),lty=1);




#########################
### Percentage Change ###
#########################
curve(pct.Fcond.rho, ylim=ylim.pct.Fcond, col=2, ylab="% of Decomposition", xlab="Quantile Index",
      main="Percentage decomposition", from=l.trim/100, to=u.trim/100);
curve(lbound.pct.cond.rho, add=T, col=2, lty=2);
curve(ubound.pct.cond.rho, add=T, col=2, lty=2);
curve(pct.Heck.Fcond.rho, add=T, col="dark red", lty=3);
curve(pct.Fcond.pi, add=T, col=3);
curve(lbound.pct.cond.pi, add=T, col=3, lty=2);
curve(ubound.pct.cond.pi, add=T, col=3, lty=2);
curve(pct.Heck.Fcond.pi, add=T, col="dark green", lty=3);
curve(pct.Fcond.beta, add=T, col=4);
curve(lbound.pct.cond.beta, add=T, col=4, lty=2);
curve(ubound.pct.cond.beta, add=T, col=4, lty=2);
curve(pct.Heck.Fcond.beta, add=T, col="dark blue", lty=3);
curve(pct.Fcond.Fz, add=T, col=6);
curve(lbound.pct.cond.Fz, add=T, col=6, lty=2);
curve(ubound.pct.cond.Fz, add=T, col=6, lty=2);
curve(pct.Heck.Fcond.Fz, add=T, col=6, lty=3);
abline(h=0);
legend("bottomright", legend=c(expression(rho),expression(pi),expression(beta),"Fz"), col=c(2:4,6),lty=1);

# Without Heckman #
curve(pct.Fcond.rho, ylim=ylim.pct.Fcond, col=2, ylab="% of Decomposition", xlab="Quantile Index",
      main="Percentage decomposition", from=l.trim/100, to=u.trim/100);
curve(lbound.pct.cond.rho, add=T, col=2, lty=2);
curve(ubound.pct.cond.rho, add=T, col=2, lty=2);
curve(pct.Fcond.pi, add=T, col=3);
curve(lbound.pct.cond.pi, add=T, col=3, lty=2);
curve(ubound.pct.cond.pi, add=T, col=3, lty=2);
curve(pct.Fcond.beta, add=T, col=4);
curve(lbound.pct.cond.beta, add=T, col=4, lty=2);
curve(ubound.pct.cond.beta, add=T, col=4, lty=2);
curve(pct.Fcond.Fz, add=T, col=6);
curve(lbound.pct.cond.Fz, add=T, col=6, lty=2);
curve(ubound.pct.cond.Fz, add=T, col=6, lty=2);
abline(h=0);
legend("bottomright", legend=c(expression(rho),expression(pi),expression(beta),"Fz"), col=c(2:4,6),lty=1);

# Percentage Change of DF #
matplot(ys[qn.plot], cbind(pct.DFcond.rho,lbound.pct.DFcond.rho,ubound.pct.DFcond.rho,pct.Heck.DFcond.rho,
                           pct.DFcond.pi,lbound.pct.DFcond.pi,ubound.pct.DFcond.pi,pct.Heck.DFcond.pi,
                           pct.DFcond.beta,lbound.pct.DFcond.beta,ubound.pct.DFcond.beta,pct.Heck.DFcond.beta,
                           pct.DFcond.Fz,lbound.pct.DFcond.Fz,ubound.pct.DFcond.Fz,pct.Heck.DFcond.Fz)[qn.plot,],
        col=c(2,2,2,"dark red",3,3,3,"dark green",4,4,4,"dark blue",6,6,6,6), type="l", lty=rep(c(1,2,2,3),4),
        ylab="% of Decomposition", xlab="Quantile",
        main="Percentage decomposition");
abline(h=0);
legend("bottomright", legend=c(expression(rho),expression(pi),expression(beta),"Fz"), col=c(2:4,6),lty=1);

# Percentage Change of DF without Heckman #
matplot(ys[qn.plot], cbind(pct.DFcond.rho,lbound.pct.DFcond.rho,ubound.pct.DFcond.rho,
                           pct.DFcond.pi,lbound.pct.DFcond.pi,ubound.pct.DFcond.pi,
                           pct.DFcond.beta,lbound.pct.DFcond.beta,ubound.pct.DFcond.beta,
                           pct.DFcond.Fz,lbound.pct.DFcond.Fz,ubound.pct.DFcond.Fz)[qn.plot,],
        col=c(2,2,2,3,3,3,4,4,4,6,6,6), type="l", lty=rep(c(1,2,2),4),
        ylab="% of Decomposition", xlab="Quantile",
        main="Percentage decomposition");
abline(h=0);
legend("bottomright", legend=c(expression(rho),expression(pi),expression(beta),"Fz"), col=c(2:4,6),lty=1);


# Combine the selection decompositions #
curve(pct.Fcond.selection, ylim=ylim.pct.Fcond, col=2, ylab="% of Decomposition", xlab="Quantile Index",
      main="Percentage Decomposition of Observed Wage", from=l.trim/100, to=u.trim/100);
curve(lbound.pct.cond.selection.2, add=T, col=2, lty=2);
curve(ubound.pct.cond.selection.2, add=T, col=2, lty=2);
curve(pct.Heck.Fcond.selection, add=T, col="dark red", lty=3);
curve(pct.Fcond.beta, add=T, col=4);
curve(lbound.pct.cond.beta.2, add=T, col=4, lty=2);
curve(ubound.pct.cond.beta.2, add=T, col=4, lty=2);
curve(pct.Heck.Fcond.beta, add=T, col="dark blue", lty=3);
curve(pct.Fcond.Fz, add=T, col=6);
curve(lbound.pct.cond.Fz.2, add=T, col=6, lty=2);
curve(ubound.pct.cond.Fz.2, add=T, col=6, lty=2);
curve(pct.Heck.Fcond.Fz, add=T, col=6, lty=3);
abline(h=0);
legend("bottomright", legend=c("selection",as.expression(bquote(beta)),"Fz"), col=c(2,4,6),lty=1);

# Combine the selection decompositions without Heckman #
curve(pct.Fcond.selection, ylim=ylim.pct.Fcond, col=2, ylab="% of Decomposition", xlab="Quantile Index",
      main="Combined Percentage Decomposition of Conditional QF", from=l.trim/100,to=u.trim/100);
curve(lbound.pct.cond.selection.2, add=T, col=2, lty=2);
curve(ubound.pct.cond.selection.2, add=T, col=2, lty=2);
curve(pct.Fcond.beta, add=T, col=4);
curve(lbound.pct.cond.beta.2, add=T, col=4, lty=2);
curve(ubound.pct.cond.beta.2, add=T, col=4, lty=2);
curve(pct.Fcond.Fz, add=T, col=6);
curve(lbound.pct.cond.Fz.2, add=T, col=6, lty=2);
curve(ubound.pct.cond.Fz.2, add=T, col=6, lty=2);
abline(h=0);
legend("bottomright", legend=c("selection",as.expression(bquote(beta)),"Fz"), col=c(2,4,6),lty=1);

# Percentage Change of DF without Heckman #
matplot(ys[qn.plot], cbind(pct.DFcond.selection,lbound.pct.DFcond.selection.2,ubound.pct.DFcond.selection.2,
                           pct.DFcond.beta,lbound.pct.DFcond.beta,ubound.pct.DFcond.beta,
                           pct.DFcond.Fz,lbound.pct.DFcond.Fz,ubound.pct.DFcond.Fz)[qn.plot,],
        col=c(2,2,2,4,4,4,6,6,6), type="l", lty=rep(c(1,2,2),3),
        ylab="% of Decomposition", xlab="Quantile", ylim=ylim.pct.Fcond, 
        main="Combined Percentage Decomposition of Conditional DF");
abline(h=0);
legend("bottomright", legend=c("selection",as.expression(bquote(beta)),"Fz"), col=c(2,4,6),lty=1);




# Quantile Function Change #
curve(QF2.cond.c1rho(x) - QF2.cond(x), ylim=ylim.QF.cond.diff, col=2, 
      ylab="Component Decomposition", xlab="Quantile Index",
      main="Decomposition of Conditional QF", from=l.trim/100, to=u.trim/100);
curve(ubound.QF.cond.joint$F2.cond.c1rho(x) - lbound.QF.cond.joint$F2.cond(x), add=T, col=2, lty=2);
curve(lbound.QF.cond.joint$F2.cond.c1rho(x) - ubound.QF.cond.joint$F2.cond(x), add=T, col=2, lty=2);
curve(QF2.Heck.cond.c1rho(x) - QF2.Heck.cond(x), add=T, col="dark red", lty=3);
curve(QF2.cond.c1rhopi(x) - QF2.cond.c1rho(x), add=T, col=3);
curve(ubound.QF.cond.joint$F2.cond.c1rhopi(x) - lbound.QF.cond.joint$F2.cond.c1rho(x), add=T, col=3, lty=2);
curve(lbound.QF.cond.joint$F2.cond.c1rhopi(x) - ubound.QF.cond.joint$F2.cond.c1rho(x), add=T, col=3, lty=2);
curve(QF2.Heck.cond.c1rhopi(x) - QF2.Heck.cond.c1rho(x), add=T, col="dark green", lty=3);
curve(QF2.cond.c1rhopibeta(x) - QF2.cond.c1rhopi(x), add=T, col=4);
curve(ubound.QF.cond.joint$F2.cond.c1rhopibeta(x) - lbound.QF.cond.joint$F2.cond.c1rhopi(x), add=T, col=4, lty=2);
curve(lbound.QF.cond.joint$F2.cond.c1rhopibeta(x) - ubound.QF.cond.joint$F2.cond.c1rhopi(x), add=T, col=4, lty=2);
curve(QF2.Heck.cond.c1rhopibeta(x) - QF2.Heck.cond.c1rhopi(x), add=T, col="dark blue", lty=3);
curve(QF1.cond(x) - QF2.cond.c1rhopibeta(x), add=T, col=6);
curve(ubound.QF.cond.joint$F1.cond(x) - lbound.QF.cond.joint$F2.cond.c1rhopibeta(x), add=T, col=6, lty=2);
curve(lbound.QF.cond.joint$F1.cond(x) - ubound.QF.cond.joint$F2.cond.c1rhopibeta(x), add=T, col=6, lty=2);
curve(QF1.Heck.cond(x) - QF2.Heck.cond.c1rhopibeta(x), add=T, col=6, lty=3);
curve(QF1.cond(x) - QF2.cond(x), add=T, col=1);
curve(ubound.QF.cond.joint$F1.cond(x) - lbound.QF.cond.joint$F2.cond(x), add=T, col=1, lty=2);
curve(lbound.QF.cond.joint$F1.cond(x) - ubound.QF.cond.joint$F2.cond(x), add=T, col=1, lty=2);
curve(QF1.Heck.cond(x) - QF2.Heck.cond(x), add=T, col=1, lty=3);
abline(h=0);

legend("bottomright", legend=c(expression(rho),expression(pi),expression(beta),"Fz","Total"), col=c(2:4,6,1),lty=1);



# Percentage Change with Joint CB #
curve(pct.Fcond.rho, ylim=ylim.pct.Fcond, col=2, ylab="% of Decomposition", xlab="Quantile Index",
      main="Percentage Decomposition of Conditional QF with Joint CB", from=l.trim/100, to=u.trim/100);
curve((ubound.QF.cond.joint$F2.cond.c1rho(x) - lbound.QF.cond.joint$F2.cond(x))/
        (lbound.QF.cond.joint$F1.cond(x) - ubound.QF.cond.joint$F2.cond(x))*100, add=T, col=2, lty=2);
curve((lbound.QF.cond.joint$F2.cond.c1rho(x) - ubound.QF.cond.joint$F2.cond(x))/
        (ubound.QF.cond.joint$F1.cond(x) - lbound.QF.cond.joint$F2.cond(x))*100, add=T, col=2, lty=2);
curve(pct.Heck.Fcond.rho, add=T, col="dark red", lty=3);
curve(pct.Fcond.pi, add=T, col=3);
curve((ubound.QF.cond.joint$F2.cond.c1rhopi(x) - lbound.QF.cond.joint$F2.cond.c1rho(x))/
        (lbound.QF.cond.joint$F1.cond(x) - ubound.QF.cond.joint$F2.cond(x))*100, add=T, col=3, lty=2);
curve((lbound.QF.cond.joint$F2.cond.c1rhopi(x) - ubound.QF.cond.joint$F2.cond.c1rho(x))/
        (ubound.QF.cond.joint$F1.cond(x) - lbound.QF.cond.joint$F2.cond(x))*100, add=T, col=3, lty=2);
curve(pct.Heck.Fcond.pi, add=T, col="dark green", lty=3);
curve(pct.Fcond.beta, add=T, col=4);
curve((ubound.QF.cond.joint$F2.cond.c1rhopibeta(x) - lbound.QF.cond.joint$F2.cond.c1rhopi(x))/
        (lbound.QF.cond.joint$F1.cond(x) - ubound.QF.cond.joint$F2.cond(x))*100, add=T, col=4, lty=2);
curve((lbound.QF.cond.joint$F2.cond.c1rhopibeta(x) - ubound.QF.cond.joint$F2.cond.c1rhopi(x))/
        (ubound.QF.cond.joint$F1.cond(x) - lbound.QF.cond.joint$F2.cond(x))*100, add=T, col=4, lty=2);
curve(pct.Heck.Fcond.beta, add=T, col="dark blue", lty=3);
curve(pct.Fcond.Fz, add=T, col=6);
curve((ubound.QF.cond.joint$F1.cond(x) - lbound.QF.cond.joint$F2.cond.c1rhopibeta(x))/
        (lbound.QF.cond.joint$F1.cond(x) - ubound.QF.cond.joint$F2.cond(x))*100, add=T, col=6, lty=2);
curve((lbound.QF.cond.joint$F1.cond(x) - ubound.QF.cond.joint$F2.cond.c1rhopibeta(x))/
        (ubound.QF.cond.joint$F1.cond(x) - lbound.QF.cond.joint$F2.cond(x))*100, add=T, col=6, lty=2);
curve(pct.Heck.Fcond.Fz, add=T, col=6, lty=3);
abline(h=0);

legend("bottomright", legend=c(expression(rho),expression(pi),expression(beta),"Fz"), col=c(2:4,6),lty=1);

dev.off()